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ABSTRACT 

We study the topology of reionization using accurate three-dimensional radiative 
transfer calculations post-processed on outputs from cosmological hydrodynamic sim- 
ulations. In our simulations, reionization begins in overdense regions and then "leaks" 
directly into voids, with filaments rcionizing last owing to their combination of high 
recombination rate and low emissivity. This result depends on the uniquely-biased 
emissivity field predicted by our prescriptions for star formation and feedback, which 
have previously been shown to account for a wide array of measurements of the post- 
reionization Universe. It is qualitatively robust to our choice of simulation volume, 
ionizing escape fraction, and spatial resolution (in fact it grows stronger at higher 
spatial resolution) even though the exact overlap redshift is sensitive to each of these. 
However, it weakens slightly as the escape fraction is increased owing to the reduced 
density contrast at higher redshift. We also explore whether our results are sensitive 
to commonly-employed approximations such as using optically-thin Eddington tensors 
or substantially altering the speed of light. Such approximations do not qualitatively 
change the topology of reionization. However, they can systematically shift the over- 
lap redshift by up to Az ~ 0.5, indicating that accurate radiative transfer is essential 
for computing reionization. Our model cannot simultaneously reproduce the observed 
optical depth to Thomson scattering and ionization rate per hydrogen atom at z = 6, 
which could owe to numerical effects and/or missing early sources of ionization. 

Key words: radiative transfer — cosmology: theory - early Universe — diffuse 
radiation — intergalactic medium 



1 INTRODUCTION 

In the concordance ACDM cosmology, the sources that dom- 
inate cosmological reionization form predominantly in over- 
dense regions. In the presence of an inhomogeneous inter- 
galactic medium (IGM), the way in which ionization fronts 
(I-fronts) propagate into the IGM depends on the spa- 
tial distributions of the gas density and sources. During 
the early stages of reionization, I-fronts proceed preferen- 
tially from the overdense knots, where sources form, to- 
ward underdense regions, where the hydrogen recombina- 
tion rates are lower; this is referred to as "inside-out" (IO) 
reionization. By contrast, in the final stages of reioniza- 
tion or within regions that have already reionized, reion- 
ization proceeds predominantly from voids towards over- 
densities; this is referred to as "outside-in" (OI) reioniza- 
tion. The way in which reionization proceeds from the ini- 
tial IO phase to the final OI phase affects the evolution of 
the size spectrum of ionized regions, leaving observable im- 
prints on th e power spectrum of fluctuations in the 21cm 
background jFurlanetto fc Ohl 120051 ; iMcQuinn etaH 12007m 
and t he galaxy-21cm cross-correlation function I Lidz et ahl 
2009). Additionally, the topology of reionization directly de- 
termines the dependence of the volume-averaged hydrogen 



recombination rate on the neutral hydrogen fraction through 
the clumping factor because it determines the "order" in 
which regions with different recombination rates are reion- 
ized. Hence, it is an important ingredient in observational 
constraints on the strength of the ionizing background as 
well as in semi-analytic models of reionization. 

For these reasons, the topology of reionization has been 
the subject of numerous recent investigations. In one of the 
pioneering radiat ive hydrodynamic simulati ons of cosmolog- 
ical reionization, IGnedin fc Ostrikerl (|l997h found that the 
clumping factor of ionized hydrogen increases monotonically 
with time. While they did not specifically discuss the topol- 
ogy of reionization, this trend suggests a completely OI 
topology because the baryonic clumping factor is smaller 
in voids than in overdensities. The approximations in this 
work effectively smoothed the ionizing background over a 
large cosmological volume, and the resulting topology was 
probably an artifact of this treatment. Similar results have 
been found in simulations that i ntroduce the ionizing back- 
ground as a boundary condition dNakamoto et al.ll200ll ) . In- 
deed, soon afterwards, IGnedin! i|2000l ) found that a more ac- 
curate treatment for the spatial distribution of the ionizing 
sources yields a more IO-like reionization in the sense that 
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the clumping factor of ionized hydrogen decreases monoton- 
ically in time. 



In the same year, iMiralda-Escude et all (|2000h formu- 
lated an influential analytic treatment for the final stages of 
reionization based on the idea that the last regions to reion- 
ize are those that combine high density with low emissivity. 
This model also describes the evolution of the ionization 
field within HII regions. However, it does not address the 
early stages of reionization, which are currently being tack- 
led by large-scale radiative transfer simulations that attempt 
to model the expansion of I-fronts out of small haloes. These 
calculations now consistently predict an IO topology in the 
sense that the ratio of the mass-averaged to the volume- 
averaged ionized hydrogen fraction xm/xv is greater than 
unity at all tim es (|lliev et al.ll2006al . l2007l ; lTrac fe Cenll2007l; 
iLee et al"1l2008h . Finally, a recent semi- numerical work spec- 
ulated that treating the spatial distribution of the hydro- 
gen recombination rate accurately could lead to a hybrid 
scenario in which overdense region s ionize first, then under - 
dense regions, and then filaments (|Choudhurv et all feoOO). 
In summary, there is as yet no consensus on how the topol- 
ogy of reionization evolves from the Dark Ages until overlap. 

We have recently developed an accurate technique for 
computing cosmological reionization that uses a moment 
method t o solve the fully time -dependent radiative transfer 
equation l|Finlator et al.l 2009). Here, we use this method 
to investigate the topology of reionization. This work is 
complementary to previous studies in that we derive the 
baryon density and emissivity fields from hydrodynamic 
simulations whose baryonic physics treatments have been 
shown to reproduce a wide range of observations of the post- 
reionization Universe , including the high-redshift galaxy 
luminosity func tion (IDave et al.| I2006T) and th e metallici- 
tiesof galaxies dFinlator fc Davg|2008h. groups (IDave et al.l 
20081). and the IGM (lOppenheimer fc Pavel 120061 . |200S| ; 
Oppenheimer et al.l l2009). Of course, we are making an as- 
sumption that the same emission sources and processes dom- 
inate during reionization; in effect, we will test this assump- 
tion by comparing to available observations. Additionally, 
the present work represents an i mprovement ove r the p re- 
liminary application presented in iFinlator et al.l l|2009h in 
three important respects: (1) we now evolve the emissiv- 
ity and density fields in time rather than assuming them to 
be static; (2) the present simulations incorporate eight times 
the cosmological volume while the underlying star formation 
rates are computed using the same mass resolution; and (3) 
these calculations begin at z — 14 rather than z = 9, so they 
account more accurately for the impact of ionizing photons 
that were emitted well before the hydrogen recombination 
time at the mean IGM density exceeded the Hubble time. 

In Section [21 we describe the cosmological simulation 
that we use as a basis for our post-processing radiative trans- 
fer calculations, and review our radiative transfer technique. 
In Section O we study the topology of reionization in our 
calculations. In Section [4] we compare our results with con- 
straints on the integrated optical depth to Thomson scatter- 
ing and the volume-averaged ionizing emissivity at z ~ 6. 
We discuss our results in Sections [5] Finally, we summarize 
our findings in Section [6] 



2 SIMULATIONS 

2.1 Cosmological Simulation 

We ran the cosmological hydrodynamic simulation that 
serves as an input to our post-processing radiative 
transfer calculation using our custom version of the 
pa rallel cosmological galaxy formation code Gadget- 
2 l|Springel fc Hernquisd[2~002r ). This code uses an entropy- 
conservative formulation of smoothed particle hydrodynam- 
ics (SPH) along with a tree-particle-mesh algorithm for han- 
dling gravity. It accounts for photoionization heating start- 
ing at z = 9 via a spatially u niform photoionizing back- 
ground (jHaardt fc M adau 200l|). Gas particles undergo ra- 
diative cooling under the assumption of ionization equilib- 
rium, where we account for metal-line cooling using the colli- 
sional ionization equilibrium tables of ISutherland fc Dopital 
(1993). Stars are formed from dense gas via a subresolution 
multi-phase model that tracks condens ation and evapora- 
tion in the interstellar medium following iMcKee fc Ostrikerl 
(1977). The model is tuned via a single parameter, t he star 
formation ti mescale, to reproduce the [ Kennicuttl (|l998al ) 
relation; see ISpringel fc Hernquistl |2003a|)" for details. Self- 
enrichment of star-forming gas and delayed feedback from 
old star particles are also treated. We account for galactic- 
scale superwind feedback using our momentum-driven out- 
flows with a normalization ao — 150 kms -1 . For fur- 
the r details on the physi c s treatmen ts in the simulations, 
see lOppenheimer fc Pavel (|2006l .l 2008). Our simulation sub- 
tends a cubical volume 16/i _1 Mpc long on each side and uses 
512 3 dark matter and star particles. It assumes a cosmology 
where fl M = 0.25, D. A = 0.75, H = 70 kms" 1 Mpc" 1 , 
cr 8 = 0.83, and Q b = 0.044. 

We now briefly discuss how our mass resolution com- 
pares to the mass scales that dominate reionization. Haloes 
that are less massive than the threshold for atomic cool- 
ing (8.3 x 10 7 (7/(1 + 0)) 1 ' 6 M S ) do not contribute sig- 
nifica ntly to reionizatio n owing to inefficient star forma- 
tion (|Wise fc C en 2009) and the early passage of Lyman - 
Werner photons (lHaiman et al.l Il997l ; lAhn et all 120091 '). 
hence we are interested only in more massive haloes. Ideally, 
we would like to resolve star formation in haloes down to 
roughly 2 x 1O 7 M0, the atomic cooling threshold at z = 14. 
In order to understand how our ionizing emissivity field re- 
lates to the underlying halo population, we have identified 
the dark matter haloes in our cosmological volume at the 
representative redshifts z — 9 and z — 6 using a spherical 
overdensity algorithm that determines the smallest radius 
out to which a halo's enclosed density falls below the virial 
density. 

We show in Figure [1] the predicted relationship be- 
tween halo mass M and star formation rate M„. For haloes 
above 1O 9 M0, both relationships follow a trend M, oc 
M 1 ' 3 . The superlinear scaling results from our momentum- 
driven outflows because the outflow mass loading factor 
r/w scales inversely with the velocity dispersion, preferen- 
tially suppressing the star format ion rates of low-mass sys- 
tems (|Oppenheimer fc Da ve 2006). At both redshifts, there 
is a turnover at low masses. Defining the turnover mass 
as the highest mass where the median M* is below 50% 
of the linearly extrapolated high-mass trend (extrapolated 
from M > 10 9 M Q ), we find that the turnover mass rises 



from 10 8 ' 3 M Q at 



9 to 10 8 ' 7 M o at 



At z > 9, 
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Figure 1. The simulated relationship between total halo mass and instantaneous star formation rate at z = 9 (upper, blue locus) and 
2 = 6 (lower, red locus). The simulated relationships roughly follow SFRoc M 13 at high masses, with normalizations and low-mass 
cutoffs that evolve in time. The red dot-dashed line indicates the current observational limit at z = 6 (see text), which translates to a 
halo mass of 10 9 ' 5 M Q . 




8 9 10 11 

log(M) 



Figure 2. (Top) The simulated halo mass functions at z = 9 (solid) and 2 = 6 (dashed). (Bottom) The corresponding ionizing luminosity- 
weighted halo mass functions in units of ionizing photons per hydrogen atom per Hubble time per mass bin. The open squares indicate 
the lower mass limits above which (right to left) 0.2, 0.4, 0.6, and 0.8 of the total ionizing luminosity is emitted. The vertical short dashed 
line denotes the 64 particle resolution limit for baryonic physics, indicating that our simulation resolves star formation in haloes more 
massiv e than 2 X 10 s Mp,. T he vertical long-dashed line indicate the halo mass where gas infall is expected to be suppressed by 50% at 
z = 6 (Okamoto et al. 2008]). The vertical dot-dashed segment in the lower panel indicates the current observational limit at z = 6 (see 
text). 
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the turnover owes to resolution effects, while at z — 6 
it owes entirely to the suppression of cooling into low- 
mass haloes (e.g.. rThoul fe Weinberg| [l996: Oka moto et al.| 



2008) by the nascent ionizing background (jHaardt fc Madaul 



20011 '). During the same interval, the normalization of the 
M — M* trend drops by roughly 0.25 dex owing to cosmo- 
logical expansion. 

The red dot-dashed segment in Figure [T] indicates the 
current observatio nal limit at z = 6 corresponding to 
-17.5 l|Bouwens et all 120071 ) . We obtained this 
limit using the relationship between star formation rate and 
rest-frame 1350 A luminosity that arises in our hydrody- 
namic simulation at z = fQ. 

The relationship in Figure [1] is a testable prediction of 
our model because it is related to the rest-frame UV lumi- 
nosity function of Lyman-break galaxies at z = 6. It has pre- 
viously been shown (albeit with slightly different parameter 
choices) th at our model nicel y reproduces th e observed nor- 
maliz ation l|Dave et al.ll2006h and evolution (|Bouwens et alj 
120071 ) of the high redshift galaxy luminosity function. Hence 
the trends in Figure [1] represent at least a plausible model 
for the sources of ionizing photons at z > 6. 

The top panel of Figure [2] shows the simulated halo 
mass functions at the same redshifts. The mass functions 
rise smoothly to lower masses until roughly 6.0 x 10 7 Mq, 
which is the 20 dark matter particle resolution limit for 
haloes and lies close to the above target resolution. Convert- 
ing the dark matter haloes in the top panel into an ionizing 
emissivity field involves a host of additional assumptions re- 
garding gas cooling and star formation. The novelty of our 
method is that these processes are already treated in great 
detail by our hydrodynamic simulation, so that we are not 
obliged to make any additional assumptions when generat- 
ing the emissivity field snapshots that we use for comput- 
ing reionization save for the choice of ionizing escape frac- 
tio n / esc . By convolving each halo's stellar population with 
the lBruzual fc Chariot! (120031 ) population synthesis models 
assuming a IChabrierl <|2003h IMF, we have computed the 
corresponding ionizing luminosity-weighted halo mass func- 
tions. We show these in the bottom panel in units of ionizing 
photons emitted per hydrogen atom per Hubble time per 
mass bin. These curves show how haloes of different masses 
contribute to reionization, and the volume- averaged ionizing 
emissivity per hydrogen atom into the IGM is simply their 
integral multiplied by the escape fraction. 

At both redshifts, the ionizing emissivity per mass bin 
increases with decreasing mass before turning over below 
10 9 A/q. The trend at high masses is natter than the mass 
function (note that the y-axes in the bottom and top panels 
span the same dynamic range) because more massive haloes 
have higher star formation rates (Figure [TJ, and it increases 
with decreasing mass because the high abundance of low- 
mass haloes wins over their low individual star formation 
rates. 

The low-mass turnover at each redshift reflects the be- 
havior in Figure [I] At z = 9, the turnover lies below 



1 Convolving our simulated stellar populations' star formation 
histories and metallicities with the Bruzual fc Chariotl d2003T ) 
models and assuming no dust, this is given by Afuv AB = 
-21.10 - 2.881og(M»), where M, is in M yr" 1 . 



2 x 10 Mq and owes entirely to low resolution since the 
thre shold for resolved st ar formation rates is 64 star parti- 
cles dFinlator et al.ll2007T ) . We crudely suggest this threshold 
with a vertical short dashed line at 64 dark matter par- 
ticle masses = 10 829 M©. By z = 6, the uniform ionizing 
back ground that we employ in the hydrodynamic simula- 
tion (Haardt & MadaukOOll ) boosts the minimum halo mass 
for efficient gas infall to ~ fewxl0 8 Mp) (|Thoul fc Weinberg] 
1 19961 ; lokamoto et al.ll2008l) . which we indicate with a verti- 
cal long dashed line. This is well above the 64-particle min- 
imum for resolved star formation histories, hence the emis- 
sivity field during the latter stages of reionization is very 
well-resolved. Note that the suppression of star formation 
in low-mass haloes at late ti mes bears some r esemblance to 
self-regulated scenarios (e.g- Jlliev et~aill2007l ) although our 
prescription for regulating star formation is different. 

The red vertical dot-dashed segment translates the ob- 
servational limit at z = 6 from Figure Q] into halo mass, and 
the squares indicate the minimum halo mass above which 
(from right to left) 20, 40, 60, and 80% of the ionizing pho- 
tons are produced. Comparing the dot-dashed line and the 
boxes suggests that, if / osc is constant, then current obser- 
vations are sensitive to only 50% of the total ionizing lumi- 
nosity density at z = 6. 

2.2 Radiative Transfer Calculations 

We extract snapshots from this simulation at roughly 40- 
Myr intervals and map their baryonic density and emissiv- 
ity fields onto a Cartesian grid for the radiative transfer 
integration. We divide the total mass associated with SPH 
particles that lie near cell boundaries between the cells by 
summing incomplete gamma functions to their equivalent 
Plummer SPH smoothing kernels. The gas temperature in 
each cell is the mass- weighted average of the temperature of 
the gas particles that lie within the cell. We omit from the 
grid those SPH particles whose density exceeds the thresh- 
old for star formation because their effect on the ionizing 
photons is implicitly accounted for in the choice of ionizing 
escape fraction. We assume that all gas is completely neutral 
at z = 14. We use radiative grids incorporating 32 3 , 48 3 and 
64 3 cells, yielding radiative spatial resolutions Axr of 500, 
333, and 250 comoving ft _1 kpc, respectively. Poisson noise 
in the derived density fields is insignificant because, even at 
our highest resolution, there are on average (512/64) 3 = 512 
gas particles within each cell. 

We compute each cell 's emissivity by convolving its stel- 
lar populations with the iBruzual fc Chariot! (|2003l ) stellar 
population synthesis models, interpolating to the correct age 
and metallicity for each star particle. We tune the uniform 
ionizing escape fraction to 13% so that the volume-averaged 
neutral fraction falls to roughly 10~ 3 at z — 6. We group all 
ionizing photons into a single frequency group at the hydro- 
gen ionizing threshold. 

During the radiative transfer computation, we update 
the total baryon densities, emissivities, and temperatures 
using new snapshots from the cosmological simulation while 
holding the radiative variables and ionization fractions con- 
stant. The radiative transfer calculation accounts for cosmo- 
logical expansion using the same cosmology as the hydrody- 
namical simulation. 

Our radiative transfer calculations solve the moments 
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of the radiative transfer equation in cosmological comov- 
ing coordinates. After each 1-Myr timestep, we suspend the 
time- dependent integration of the moment equations and 
use a (time-independent) ray-tracing technique to compute 
the full angular dependence of the radiation field. From 
this, we derive the updated Eddington tensor field, which 
is needed to close the moment hierarchy. We use a number 
of optimizations to render this t echnique comp u tation ally 
feasible, as described in detail in iFinlator et all l|2009h . In 
brief, (i) we only update a computational cell's Eddington 
tensor when its photon number density has changed by more 
than fj = 5% since its Eddington tensor was last updated; 
(ii) we terminate a ray-tracing computation when the opti- 
cal depth between the source and the target cell along the 
ray exceeds 6; and (iii) we switch from using one to using 
two layers of periodic replica volumes in order to mimic the 
effects of periodic boundaries once the volume-averaged neu- 
tral fraction drops below 50%. In IFinlator et all (|2009l ), we 
used an extensive suite of parameter convergence tests to 
verify that each of these optimizations introduces errors in 
the ionization fractions of at most 10%. 

In this work, we introduce an additional optimization 
that enables the code to transition smoothly between the 
different ray-tracing approaches that are appropriate for the 
optically thick and thin regimes. Before cosmological H II re- 
gions overlap, the highly inhomogeneous opacity and emis- 
sivity fields lead to strong spatial variations in the Eddington 
tensor field. These can only be treated accurately by using 
ray-tracing to compute the optical depths from every posi- 
tion to every source. However, as reionization proceeds, a 
growing fraction of the computational cells lie deep within 
large HII regions where the Eddington tensors are dom- 
inated by other sources located within the same HII re- 
gion and towards which the optical depth is negligible. In 
these regions, simply assuming that the optical depth to 
every source is zero becomes increasingly valid, and the 
computationally-efficient optically thin approximation is ap- 
propriate. Optimization thus involves defining a criterion 
to determine when a computational cell lies within a large 
HII region. We do this as follows: At all times, we store the 
optically thin Eddington tensor field (that is, the result from 
assuming that all optical depths are zero) in memory. Wher- 
ever the (time-dependent) photon number density exceeds 
the photon number density from the (time-independent) op- 
tically thin approximation by a factor of 2, we use the op- 
tically thin Eddington tensors rather than performing ray- 
tracing. We have used a paramet er convergence test similar 
to those in IFinlator et alj |2009l ) to verify that this choice 
incurs no more than 10% errors in the ionization fractions 
while speeding up the computation by roughly a factor of 2. 

Through extensive testing, we have found that relax- 
ing our accuracy criteria so that we recompute the Edding- 
ton tensors whenever the photon number density changes by 
10% and switch to the optically thin approximation once the 
time-dependent photon density exceeds the optically thin 
value by 10% introduces negligible changes into the topol- 
ogy and redshift of reionization (Figure [9]) although it does 
introduce typical errors of 20% into the ionization states of 
individual cells. We refer to this as the "fast" scheme and 
employ it for tests of systematics. 

We evolve the ionization field using either implicit or 
explicit finite-differencing techniques depending on the local 



ionization timescale. We do not evolve the gas temperature 
because we do not solve for the hydrodynamic response of 
the gas to the passage of ionization fronts. At each timestep, 
we iterate between the updates to the ionization and radia- 
tion fields until they converge to 10 -4 . 

Our highest-resolution simulation required roughly 
10,000 CPU hours on a shared-memory machine with 8 2.5- 
GHz Intel Xeon CPUs. The overall computation time for 
an individual reionization simulation scales with the total 
number of cells N ce n s approximately as -/V* cl 5 ls , reflecting the 
fact that th e time for the Eddin gton tensor updates scales 
in this way (IFinlator et alj|2009h . 



3 THE EARLY REIONIZATION OF VOIDS 

3.1 The Volume- Averaged Neutral Fraction 

In the top panel of Figure [3] we show how the volume- 
averaged neutral fraction evolves in time at our highest spa- 
tial resolution (Axr = 250/i~ 1 kpc). The solid curve gives 
the average over our entire simulation volume while the dot- 
ted, short-dashed, and long-dashed curves correspond to un- 
derdense, overdense, and mean-density regions, respectively. 
Note that using different ranges does not change the qual- 
itative results. Hereafter, we refer to underdense regions as 
"voids" and mean-density regions as "filaments". 

The first thing to notice is that overlap (defined as 
the point where the volume-averaged neutral fraction ihi.v 
drops below 10 -3 ) occurs by z — 6. This is a nontriv- 
ial result, given that our simulations have previously been 
shown to reproduce a wide array of observations of the post- 
reionization Universe (Section 12. 2[) and that we have not in- 
troduced any new baryonic physics in order to bring about 
reionization other than the (reasonable) choice of ionizing es- 
cape fraction. This supports previous suggestions that reion- 
ization could have been dominated b y ordinary (that is, 
not P opulation III) star formation fe.g.. lShull fc VenkatesarJ 

boost) . 

Examining the reionization history in different density 
bins, we find that overdense regions reionize first because 
they host the bulk of the ionizing sources. Around z = 9, 
photons bypass the filaments and flow into the voids, which 
ionize rapidly owing to their low recombination rates. Fil- 
aments ionize last because they possess fewer sources than 
overdense regions but higher recombination rates than voids. 
We refer to this topology, in which the voids are ionized 
significantly before the filaments, as the "inside-outside- 
middle" topology (IOM). Following overlap, the neutral hy- 
drogen fraction increases with overdensity as expected in the 
presence of a uniform ionizing background. 

In the bottom panel of Figure [3] we show how the 
volume-averaged mean ionizing intensity evolves in the same 
regions. Overdense regions see the strongest ionizing back- 
ground owing to their proximity to sources. The mean inten- 
sity in filaments always represents an average of the inten- 
sities in reionized regions close to sources and self-shielded 
regions farther away (Figure |5), but broadly it lies between 
the intensities in overdense and underdense regions. The in- 
tensity in voids is negligible until the first I-fronts bypass 
filaments at z = 9. Afterwards, voids ionize rapidly owing 
to their low recombination rates and their mean intensity 
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Figure 3. (Top) The volume-averaged neutral fraction in the IGM for the Axr = 250/i — 1 kpc calculation as a function of the age of 
the Universe (bottom axis) and redshift (top axis) overall and for three bins in overdensity. (Bottom) The mean intensity of the ionizing 
background at 912 A for the same density bins. The "notches" at roughly 40-Myr intervals (especially clear in the bottom panel) indicate 
the ti mesteps where the ioniza tion and density fields are updated during the integration. The triangle indicates the observed upper 
limit l|Bolton fc Haehneltll2007h . Voids reionize before mean-density regions owing to the spatial distribution of recombination rates and 
ionizing intensities. 



rapidly reaches the volume average. Following overlap at 
z ~ 6, the ionizing background is very nearly homogeneous. 

In Figure 31 we show how the median redshift of reion- 
ization varies with overdensity and Axn. We constructed 
this figure by determining the first redshift at which each 
cell's neutral hydrogen fraction drops below 10~ 3 . We de- 
fine each cell's baryonic overdensity at the final timestep 
(z = 5.71). Broadly, our finding that reionization occurs 
sooner in voids than in filaments is insensitive to As;r. De- 
creasing Axr. delays reionization in overdense regions by en- 
hancing I-front trapping while accelerating the reionization 
of voids. Defining the overlap epoch as the redshift at which 
the volume- averaged neutral fraction drops below 10~ 3 , we 
find that increasing the spatial resolution from Axb. = 500 
to 250/i~ kpc causes overlap to occur earlier by Az — 0.2. 
We will argue in Section 13.31 that both of these effects owe 
to more effective "tunneling" of photons through soft spots 
in the IGM at higher spatial resolution. 

The solid and short dashed - long dashed curves in Fig- 
ure [4] show how reionization proceeds when we use Edding- 
ton tensors that are derived accurately versus in the opti- 
cally thin limit, respectively. Using optically thin Edding- 
ton tensors systematically accelerates reionization in voids 
while delaying it in filaments. Overdense regions remain un- 
affected. The artificially late reionization in filaments delays 
the overlap epoch by Az — 0.2. Note that this constitutes 
the first evaluation of the consequences of incorporating in- 
accurate Eddington tensors into calculations of cosmological 
reionization. 

In order to show more intuitively how the topology in 
Figure|4]arises, we show in Figure[5]maps of overdensity (a), 



ionizing emissivity (c,f), neutral hydrogen fraction (d,g), and 
the intensity of the ionizing background (e,h) at two repre- 
sentative redshifts as well as the redshift of reionization (b) 
in the same thin slice through our cosmological volume (see 
caption for details). We note that the simulation from which 
we extract these maps is identical to the runs in Figure [4] 
except that it uses 96 3 radiative transfer cells for a radiative 
transfer resolution of Axr = 177ft ~ kpc. Because this run 
became prohibitively slow after the neutral fraction dropped 
below xhi,v < 10%, we completed its final stages using op- 
tically thin Eddington tensors. This impacts the topology of 
reionization negligibly — in fact, we found that this simula- 
tion continues all of our trends with respect to resolution — 
but for consistency we omit it from the other figures in this 
paper. 

Panel 5(b) can be compared to Figure 1 of iTrac et al.l 

(2008), who found a purely IO reionization topology. Com- 
paring panels a, c, and f reveals that, while sources lie at 
the intersections between filaments as expected, they do not 
trace the filaments smoothly out to the mean density owing 
to the low-mass cutoff in the halo mass to light ratio (Fig- 
ure [2}. This suggests that, for the most part, the filaments 
cannot self-ionize. Panels d, e, g, and h show that some fil- 
aments do ionize rapidly owing to proximity to sources, but 
at distances of more than ~a few h~ 1 Mpc from sources, 
filaments remain self-shielded even after most of the voids 
have already reionized. The tendency for most filaments to 
reionize late gives rise to our characteristic IOM reionization 
topology. 

A close examination of the emissivity maps (panels [5(c) | 
and 5(g) I shows that our simulated emissivity varies rapidly 
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with position. This is a consequence of the fact that the ion- 
izing luminosity of a single stellar population varies by 5 or- 
ders of magnitude during its first 100 Myr. Our limited mass 
resolution leads to shot noise in the number of star particles 
per grid cell whose age lies within this range, which in turn 
gives rise to the noisy emissivity field. To gauge the impact 
of this, we tried computing and using instantaneous star for- 
mation rates directly from gas particles as the source of ion- 
izing flux, rather than the star formation history from star 
particles. This would be a more appropriate approach in the 
limit that all ionizing photons are emitted instantaneously. 
We left the escape fraction constant and tuned the ionizing 
luminosity per unit star formation rate to 2.45 x 10 53 s _1 
per Mq yr _1 in order to match the volume-averaged emissiv- 
ity resulting from the star particlefl We show the resulting 
reionization topology and overlap redshift with the dotted 
green lines in Figure [9] Comparing these trends with the 
solid black lines, which reflect our fiducial emissivity pre- 
scription, we find that the two approaches yield nearly iden- 
tical results. This indicates that our overall results are not 
sensitive to stochasticity in the star formation algorithm. 

3.2 The Ratio of Ionized Fractions 

Another way to study the topology of reionization is to con- 
sider the ratio of the mass-averaged to the volume-averaged 
ionized hydrogen fractions, xm/xv, which can be thought 
of as the aver age density of ion ized regions in units of the 
mean density (|lliev et alj2006ah . For reionization in a homo- 
geneous IGM, xm/xv = 1 at all times. In an inhomogeneous 
IGM, 01 reionization implies that xm/xv ^ 1 because ion- 
izing voids increases xhii,v more rapidly than xb_u,m- Sim- 
ilarly, IO reionization implies that xm/xv ^ 1 because ion- 
izing overdense regions increases ihii.m more rapidly than 
3;hii,v- In Figure [S] we show how xm/xv evolves in time for 
three different choices of Air in the radiation grid. For the 
period z ^ 9, we smooth the simulated results with cubic 
polynomial fits in order to suppress artifacts from low time 
resolution in our set of precomputed density and emissiv- 
ity grids. The inset panel shows the original trend for the 
Air = 250/i _1 kpc calculation. 

Broadly, im/iv evolves from > 1 at early times when 
the ionized volume fraction remains low, to < 1 when the 
voids reionize, to unity once the m ajority of the universe has 
reionized (see also llliev et al ]|2007i). The early evolution is as 
expected if reionization begins ea rlier in overdense regions 
than in the Universe as a whole (|Gnedinll2000l; [ii iev et al.l 
2006a), but the tendency for xyi/xy to drop below unity 
significantly before reionization completes has not been dis- 
cussed previously. The evolution of %/iv is sensitive to the 
choice of Air: Increasing spatial resolution increases xm/xv 
at earlier times while decreasing it at late times. Physically, 
higher spatial resolution allows the simulation to better re- 
solve the high recombination rates in the overdense regions 
where reionization begins, leading to more effective trap- 
ping of I-fronts at early times. At late times, increasing the 
resolution allows the simulation to resolve the tendency for 

2 This is 2.64 times the iKennicutj jl998ah relation. The extra 
ionizing luminosity can be attributed to the low metallicity at 

2^6. 



I-fronts to "leak out" through the porous IGM into voids. 
As the bottom panel of Figure|6]indicates, this leaking effect 
is so efficient that, at our highest spatial resolution, xm/xv 
drops below unity around the time when ihi.v drops below 
50% (at z m 7). In other words, by the time reionization is 
halfway complete, the average ionized region is underdense. 
This is a direct consequence of the late reionization of fila- 
ments. 

3.3 The Mean Free Path of Ionizing Photons 

The tendency of ionizing photons to stream preferentially 
in directions where the IGM is less dense naturally leads to 
the formation of ionized "tunnels" that connect sources with 
voids. If this tunnel formation indeed dominates the topol- 
ogy of reionization, and if the tunnels are small compared 
to Air, then the mean free path of ionizing photons Ampp 
from sources should increase with decreasing Air,. This is 
because, as the tunnelling process is better-resolved, an in- 
creasing fraction of ionizing photons travel directly from 
sources through tunnels into voids rather than being artifi- 
cially absorbed nearby. To test this idea, we have computed 
Amfp as a function of redshift by casting rays from each 
source (where a source is a computational cell with nonzero 
emissivity) in 1280 directions that uniformly sample the unit 
sphere and then computing the distance travelled until the 
optical depth exceeds 6. By repeating this exercise for each 
combination of Air and redshift, we follow the growth of 
Amfp and its dependence on Air. 

We show the resulting trends of mean free path versus 
ionized volume fraction in Figure [7] During the early stages 
of reionization, when the neutral fraction is above 99%, 
higher resolution leads to higher Amfp- This is the result 
that we expected: At higher spatial resolution, the tendency 
of photons to bore holes through the IGM from sources to 
voids is better resolved, leading to a higher Amfp at a given 
ionization state. In fact, Figure [7] indicates that even our 
highest-resolution computation does not fully capture all 
the relevant substructure, suggesting that the characteris- 
tic size of the tunnels is smaller than 250 comoving /i _1 kpc. 
It is likely that using a spatial resolution that is comparable 
to the virial radius of the dominant haloes is necessary for 
full resolution convergence. Unfortunately, this requirement 
is roughly a factor of 10 higher than what we have achieved 
here, and would require either a much smaller cosmologi- 
cal v olume (thus increasing cosmic variance and biasing ef- 
fects; iBarkana fe Loebll2004h , an adaptive RT mesh, or a less 
accurate technique. Nevertheless, the trend indicates that 
reionization proceeds rapidly from overdensities into voids 
before mopping up the filaments. It also reinforces the need 
for spatial resolutions that are much higher than 1 comoving 
/i~ 1 Mpc when studying the topology of reionization. 

As reionization proceeds, the trend of Amfp versus reso- 
lution inverts so that higher resolution leads to lower Amfp- 
This occurs as the topology changes from one in which ion- 
ized tunnels thread an otherwise opaque IGM to one in 
which increasingly isolated overdensities are separated from 
each other by regions that are already reionized. The remain- 
ing neutral regions dominate the IGM opacity and shrink 
as reionization proceeds. Higher spatial resolution prevents 
them from being averaged with the low recombination rates 
in neighboring voids. At this stage, reionization proceeds 
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Figure 4. The median redshift of reionization as a function of overdensity p/p (curves) and overlap redshift (horizontal lines) for different 
choices of A%. The short dashed - long dashed curve results from using Axr = 250/i -1 kpc along with optically thin (OT) Eddington 
tensors. Broadly, the early reionization of voids is robust to the choice of Axr as well as the accuracy of the Eddington tensors. Increasing 
the spatial resolution allows us to resolve both the high recombination rates in overdense regions and the tendency for photons to bypass 
filaments and reionize voids early. The OT approximation delays overlap by Az = 0.2 by delaying reionization in filaments. 



from the voids back into regions of increasing densit y and 
decreasing spatial scale ( Miralda-Escude et al.l l200C0. Our 
simulations suggest that this topology could dominate by 
the time the neutral fraction drops below 50%, and possibly 
earlier. 



4 COMPARISON TO OBSERVATIONAL 
CONSTRAINTS 

The goal of the present work is to take a hydrodynamical 
model that has enjoyed considerable success in accounting 
for observational constraints from the post-reionization Uni- 
verse and to ask how reionization proceeds in this model. Be- 
cause we have introduced no extra physical assumptions save 
for the choice of ionizing escape fraction and have tuned this 
to achieve reionization by z = 6, it is of interest to compare 
our results with additional observational constraints from 
the reionization epoch, specifically Thomson optical depth 
measurements from the Wtlkinson Microwave Anisotropy 
Probe (WMAP), and observations of the Lyman-Q forest. In 
this Section we show that (1) our simulated optical depth to 
Thomson scattering underproduces the observed value and 
(2) our simulations overproduce the mean ionizing intensity 
at z — 6. We discuss the implications of these findings. 

4.1 The Integrated Optical Depth to Electron 
Scattering 

In Figure [S] we show the variation of the integrated op- 
tical depth to Thomson scattering r cs with redshift. The 
thick solid curve corresponds to pure hydrogen reionization 



at our highest spatial resolution and accuracy, and indicates 
an integrated optical depth of 0.047, 0.023 below the 1-er 
observational constraints. Accounting for Helium reioniza- 
tion (thick dotted curve; see caption for details) boosts r ea 
to 0.051. The discrepancy between the observed and simu- 
lated r cs , which reflects the history of the volume- averaged 
electron number density, suggests that reionization occurs 
too suddenly or too late in our calculations. 

A number of effects could explain this discrepancy. On 
the numerical side, our hydrodynamic simulation may not 
resolve all of the relevant star formation at early times (Fig- 
ure [2J • Increased star formation at early times could start 
reionization earlier and boost r cs without changing xu i,v at 
late times, when the ionizing b ackground is dom inated by 
more massive haloes (Figure [U llliev et al.ll2007h . The role 
of early star formation could be further enhanced through 
a more detailed treatment of the ionizing escape fraction, 
/esc- Currently, we have set / osc to a time-invariant value 
of 13% in order to obtain the end of reionization at z ~ 6. 
However, this is a free parameter a nd could be larger at ear- 
lier times (e.g., IWise fc Cenl 12009). Raising / csc to 1 (thin 
curves) brings our simulated optical depth into agreement 
with the observed 68% confidence interval, but makes the 
redshift of overlap z — 7.90. Note that this change slightly 
modifies the qualitative topology of reionization, but does 
not make it purely IO (Figure [5J. Third, the finite spa- 
tial resolution affects r cs because reionization occurs more 
rapidly at higher Axr (see Section [3j|. For radiative transfer 
grid resolutions of Ax = 500, 333, 250/i~ 1 kpc, ihi.v drops 
below 10 -3 at z = 5.97, 6.05, and 6.20, respectively. Hence 
the choice of Axr, introduces an uncertainty of approxi- 
mately Az —0.2. Finally, our small cosmological volume de- 
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(f) Emissivity, z = 6.49 (xhi,v = 0.1). (g) £hi,m at z = 6.49. (h) Photon density at z = 6.49. 

Figure 5. Maps of overdensity (a), redshift of reionization (b), ionizing photon emissivity (c,f), neutral hydrogen fraction (d,g), and 
the intensity of the ionizing background (e,h) in a volume slice 16h~ 1 Mpc to a side and 250h~ x kpc thick. Figures c— h are taken at two 
representative redshifts as indicated. The knots between filaments tend to host sources, but the filaments themselves have negligible 
emissivity, with the result that they reionize last (red color in panel b). 
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Figure 6. The ratio of the mass- weighted to volume- weighted ionized hydrogen fractions im/^v a- s a function of the age of the Universe 
(lower axis) and rcdshift (top axis). The top panels shows how 1^/% drops more rapidly at higher spatial resolution owing to more 
efficient breakout of I-fronts into voids. We use polynomial fits to the simulated trends for z 9. The inset panel shows the full simulated 
result for Aier = 250/i -1 kpc, where the jumps owe to the low time resolution in our set of precomputed density and emissivity fields. 
The bottom panel expands the y-axis from the top panel about xm/^v = 1 in order to show how xm/^v drops below unity well before 
rcionization completes, with the effect growing stronger at higher resolution. 
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Figure 7. The ionizing mean free path Amfp from sources as a function of the volume-weighted ionized hydrogen fraction, computed 
using three different values of Axr. The larger panel uses 6th-order polynomial fits to the simulated trends for :ehii,v <10%, while 
the inset panel shows the full original curve for Air = 250/i — 1 kpc. At early times, higher resolution leads to higher Amfp because 
photons "tunnel" more effectively from sources to voids. At late times, higher spatial resolution leads to lower Amfp because the small 
overdensities that dominate the IGM opacity are better resolved. 
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Figure 8. The integrated optical depth to Thomson scattering t os as a function of the age of the Universe (bottom axis) and redshift 
(top axis). Thick curves correspond to Axr = 250/i kpc, / eac = 0.13, and fiducial accuracy while thin curves correspond to Ax^ = 
SOO/i^kpc, / csc = 1, and the fast scheme (Section l2,2l l, Solid curves show the result from pure hydrogeon reionization, and dotted curves 
include the contribution of Helium assuming that Helium is doubly ionized after 2 = 3 and singly ionized with riHell/ n Hc = n Hll/ n H f° r 
z > 3. The dashed lines indicate the 68% confidence intervals for r cs arising from combining WMAP-5 with distance m easurements from 
Type la supernovae and baryon acoustic oscillations, t cs — 0.084 ±0.016 llffinshaw et alj|2009l ; iKomatsu et al.i r2009). Our simulations 
do generate enough ionizing photons to match the observed constraint on r cs , but our fiducial choice of / csc may underproduce it. 



lays reionization by A 2 = 0.1 because of the lack of long- 
wavelength density fluctuations (|Barkana fc Loebll2004) and 
introduce s a random uncertainty of Az = 1 owing to cosmic 
variance (Bar kana fc Loebl 12004 llliev et all l2006al ). which 
can be translated to an effect on r cs via d,T es /dz le i on ~ 0.008. 
In short, the systematic uncertainties in our reionization his- 
tory owing to numerical effects and assumptions are suffi- 
cient to account for the missing optical depth. 

But there are fu rther observational uncerta inties related 
to the observed T es ■ IShull fc~V enkatcsan ( 20oH) have argued 
that the systematic uncertainties in computing the residual 
electron fraction leftover after recombination lead to sys- 
tematic uncertainties in r cs of order ~ 0.01. If true, then 
this could reduce the amount of Thomson scattering that 
galaxies are responsible for producing, bringing our simula- 
tions into better agreement with observations. The unknown 
time depe ndence of reionization also renders r cs uncertain 
(although iMortonson fe Hull2008l have argued that realistic 
reionization histories systematically increase r os over the in- 
stantaneous reionization value, which would further increase 
the discrepancy with our simulations). 

Finally, it is also possible that reproducing r cs 
in large cosmological volumes requires additional phys- 
ical processes such as self-regulation of star formation 
in lo w mass haloes £llicv ct al. 2007), Population III 
stars dTrac fc Cenll2007l;IShin et al.ll2008f) . X-rays from e arly 
black holes |Shull fc Venkatesanl l2008f) or supernovae l|Onl 
l200ll ) , or primordial magnetic fields (fSchleicher et al.ll2008h . 
In future work, we will use radiative hydrodynamic simula- 
tions of reionization to study self-regulation; however, for the 



present, our goal is to study how reionization follows from 
our existing baryonic physics treatments with no additional 
assumptions. 

4.2 The Hydrogen Ionization Rate at z = 6 

The solid blue curve in the bottom panel of Figure [3] 
shows how the volume- averaged hydrogen ionization rate 
per hydrogen atom T_i2 = Thi/10~ 12 varies with red- 
shift in our calculation. At z = 6, we find log(r_i2) = 
1, whereas observations of the Lyma n-q forest indicate 
log(r_i2) < —0.6 (triangle with a rrow; iBolton fc Haehnelj 
l2007l ; ISrbinovskv fc Wvithd |200ST ). Relatedly, our fiducial 
simulation yields a volume- averaged neutral fraction at 2 = 
6 of 2.4 x 10~ 6 , 50-100 t imes smaller than the observed 2- 
3xl0" 4 (|Fan et al.ll200o1 '). 

These large discrepancies could result either from / esc 
being too high or from the opacity in the reionized IGM be- 
ing too low. The fact that reionization occurs more rapidly 
at higher spatial resolution (Figure [4| indicates that we are 
forced to choose / csc artificially large in order to achieve 
overlap by 2 = 6. Meanwhile, the fact that the IGM opac- 
ity at a given ionization state increases with decreasing Axn 
(Figure^ indicates that our simulated IGM opacity at 2 = 6 
is artificially low. Both of these effects indicate that increas- 
ing our spatial resolution or incorporating a subgrid treat- 
ment for ph otoionization of structures below our spatial res- 
olution (e.g- ICiardi et al.ll2006| ; [Trac fc Cenl2007h would im- 
prove the agreement with the observed F_i2. 

We may estimate how much our F_i2 would improve 
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with resolution by considering the likely impact of structures 
below our spatial resolution. Under ionization equilibrium, 
the emissivity and ionization rate are related by 

T) = X r/47ra. (1) 

Here, jj = 2.694 x l(T 24 (/ CS c/0.13) is our simulated volume- 
averaged comoving emissivity at z = 6 in s _1 cm -3 Sr _1 ; 
X is the opacity to ionizing photons in cm -1 ; and a is the 
neutral hydrogen cross section at the Lyman limit in cm -2 . 
At Z = 6, the IGM opacity is dominated by Lyman limit 
systems, hence we may approximate \ as the reciprocal 
of the mean free path to Lyman l imit systems. Extrapo- 
lating the number density found by IStorrie-Lombardi et ail 
(1994 ), this is 84 comoving M pc in our cosmology at z = 6. 
iFaucher-Giguere et all ((2008) have estimated the additional 
impact of structures at lower column densities, finding that 
the mean free path is roughly 85[(1 + z)/4]~ 4 proper Mpc 
(or 63 comoving Mpc at z = 6). Folding these estimates 
into Equation [T] we find that resolving Lyman limit sys- 
tems could reduce T_i2 from ~ 10 to 2 . This is still larger 
than the observed upper limit o f 0.3 (|Bolton fc Haehneltl 
120071 ; ISrbinovskv fc Wvithel 120081 ). suggesting that our es- 
cape fraction may still be too large. However, given that in- 
creasing the spatial resolution also accelerates reionization 
(Figure |4|, this observation reinforces our previous conclu- 
sion that decreasing both Axr and / esc without changing 
our star formation prescription would yield improved agree- 
ment with the observations. 

Of course, while it is possible that these discrepancies 
owe to a combination of numerical effects and systematic un- 
certainties in the observed optical depth to Thomson scat- 
tering, a more interesting possibility is that they indicate a 
failing of the input physics. For example, detailed compar- 
isons with observations of galaxies and the IGM in the post- 
reionization Universe have prompted us to assume that low- 
mass galaxies generate strong outflows (Section [TJ . These 
outflows dramatically suppress the ionizing emissivity at 
early times, leading to a relatively late reionization epoch 
and a low T es . As a reference point, in our simulations, 90% 
of the ionizing photons that are emitted by z — 6 are emit- 
ted after 2 = 9. Weaker outflows at early times would boost 
early star formation rates and hence r cs . A stronger scal- 
ing between io nizing luminosity and me tallicity over what 
is found in the lBruzual fc Charlotl (|2003l ) models would en- 
hance the ionizing light to mass ratios, as would including a 
treatment for an evolving initial mass function (i.e., Popu- 
lation III star formation), or the inclusion of mini-quasars. 
Finally, an ionizing escape fraction that decreases with de- 
creasing redshift or increasing halo mass would enhance the 
ionizing luminosity into the IGM without changing the star 
formation history. However, until the numerical issues are 
more fully investigated, we cannot make robust conclusions 
about the need for new or enhanced sources of ionizing pho- 
tons at the earliest epochs. 

In summary, the failure of our simulations to reproduce 
simultaneously the observed optical depth to Thomson scat- 
tering and the hydrogen ionization rate at z — 6 could owe 
to observational uncertainties, numerical effects, or inade- 
quate physics treatments. On the other hand, given that our 
simulations have previously been shown to reproduce a wide 
variety of observations of the post-reionization Universe and 
that the only additional physical assumption that we have 



introduced is the choice of ionizing escape fraction, the level 
of agreement that we find is encouraging. Our results sup- 
port the idea that the processes governing star formation 
befo re and after reion ization may not have been very differ- 
ent jDave et al.ll2006l ). 



5 DISCUSSION 

The tendency for I-fronts to bypass filaments and escape 
directly into voids occurs generically in calculations that 
as sume static emissivi ty and density fields (as we did 
m iFinlator etafl 120091 ) . Given high enough spatial reso- 
lution, such calculations invariably pro duce H II regions 
with "butterfly w i ng" morphologie s fe.g.. lCiardi et aL I l200ll ; 
llliev et al.ll2006bl ; iTasitsiomil |2006| ) . In addition, this effect 
has recently been found in radiative hydrody namic sim- 
ulations of the first sta ges of star formation (Abe l et al.l 
120071 ; IWise fc Abel|[2008l ). This was on much smaller scales 
(< l/i _1 Mpc) than the cosmological scales that we are in- 
terested in. Nonetheless, one may understand the result in 
both cases as a competition between the timescale for pho- 
tons to bore tunnels through the IGM that connect sources 
with voids on the one hand versus the timescale for filaments 
to self-ionize and the timescale for I-fronts to burn directly 
through filaments on the other. If photons leak into voids 
more rapidly than the filaments can either be ionized or 
self-ionize, then the voids will reionize before the filaments. 

These timescales are in turn governed by three factors: 
(1) The emissivity field's bias with respect to the baryon 
density field; (2) the time evolution of the emissivity field; 
and (3) the spatial resolution. The impact of the emissiv- 
ity bias on the topology of reion ization has been explored 
elsewhere (M cQuinn et alj|2007h : for our purposes it suf- 
fices to note that a highly biased emissivity field promotes 
the late reionization of filaments by preventing them from 
self-ionizing. The time evolution of the ionizing emissivity 
matters because if the emissivity is relatively high at early 
times, then the late reionization of filaments is suppressed 
even in the presence of a biased emissivity field owing to the 
lower density contrast in the IGM. As a demonstration of 
this effect, we show the reionization topology resulting from 
setting /esc = 1 in Figure [9] (cyan dot-dashed curve). Com- 
paring this curve to the topology in our fiducial case (solid 
black curve) reveals that the higher reionization redshift 
is indeed associated with a flatter trend of median reion- 
ization redshift versus overdensity at low densities. Finally, 
high spatial resolution also promotes the early reionization 
of voids by increasing the density contrast. The topology 
presented in this paper results from a particular combina- 
tion of emissivity bias and evolution that has not been con- 
sidered previously, and it is most clearly visible only at our 
hig hest spatial resoluti on (Figure [6} • 

llliev et all l|2006al ) ran post-processed radiative trans- 
fer calculations of cosmological reionization on snapshots 
obtained from a pure N-body simulation. In contrast to 
our results, they found a purely IO topology. Given that 
their spatial resolution was comparable to ours, the differ- 
ence between our results must owe to the different emissiv- 
ity fields. In constructing their emissivity field, they con- 
sidered only haloes more massive than 2.5 x 1O 9 M0, which 
implies an even more biased emissivity field than ours (Fig- 
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ure [2} . IChoudhurv et al.l (|2009| N ) suggested that their purely 
IO topolo gy might r e sult fro m using high emissivities. In- 
deed, the llliev et al ] <|2006al ) simulation achieved reioniza- 
tion by z ~ 11, only pa 260 Myr after the first sources formed 
at z — 21, whereas our simulations are tuned to produce 
reionization at z = 6, roughly 680 Myr after the compu- 
tations begin at z — 14. We conclude that their purely IO 
topology owed to their high reionization redshift and the 
resulting lower IGM density contrast. 

Subsequent calculations involving more extended reion- 
ization histories and lower overlap reds hifts have also re- 
sulted in more n early IO-like topologies l|Trac fc Cenll20o"7l ; 
llliev et alj|2007h . hence the overlap redshift alone is not 
the on ly cause of our topology. In the case of iTrac fc Cenl 
(|2007l V we ascribe the difference to our more biased emis- 
sivity field. For example, comparing our Figure [2] with 
their Figure 1, we find that, at z = 6, the minimum halo 
mass above which 60% of ionizing photons are emitted is 
10 9,4 Mq for our simulation versus 1O 9 M0 for theirs. Their 
more extended reionization history also allows for much of 
reionization to take place at higher redshift, when density 
contrasts are lower. Our less gradual reionization history 
owes directly to our star formation prescription, which sup- 
presses the early star formation rate density with respect 
to models that do not include st rong outflows (e.g., Figure 
3 of lOppenheimer fc Pavel [200^ and leads to a more sud- 
den overlap epoch because of the strong dependence of star 
formation rate on halo mass (Figure [TJ. 

The emissivity field of llliev et ail (120071 ) is also less bi- 
ased than ours, but this time in two respects. First, they 
assume that haloes with masses between 10 8_9 Mq are the 
sites of Population III star formation, leading to lower ra- 
tios of mass to ionizing luminosity than more massive haloes. 
Second, they assume that the ionizing luminosity is propor- 
tional to halo mass, which, even in the absence of Population 
III stars, would yield a less-biased emissivity field than ours 
since our ionizing luminosity scales as A/ 1 ' 3 (Figure[T]). With 
these assumptions, it is not surprising that their filaments 
are able to self-ionize much more efficiently than ours. 

The simulations o fll iev et al. I (|2007l . l2006al ) assume that 
the speed of light is infinite whereas ours do not; could this 
assumption impact the topology of reionization? As this 
question has not been investigated in a fully numerical con- 
text previously, we have re-computed cosmological reioniza- 
tion using Axb. = 500/i _1 kpc assuming three different val- 
ues c' for the speed of light: c/10, c, and 10c. To save time, 
we perform these integrations using our less-accurate "fast" 
scheme for updating the Eddington tensors (Section 12.21) . 
Figure [9] shows how the trend of median redshift of reion- 
ization versus overdensity (curves) and the overlap redshift 
(horizontal lines) vary with the speed of light. The light solid 
curve reflects our fiducial accuracy and is copied from Fig- 
ure^ It essentially overlaps the dark solid curve, indicating 
that the impact of the somewhat less accurate Eddington 
tensors in our fast scheme is small compared to the system- 
atic offset from using an inaccurate speed of light. The me- 
dian redshift of reionization is roughly Az = 0.1 (50 Myr) 
earlier at all densities for c' = 10c and 0.2-0.4 (150 Myr) 
later for c' = c/10, although in both cases the effects are 
stronger in the underdense than in the overdense regions. 
Overlap (defined as the redshift when xhi,v drops below 
10" 3 ) occurs Az = (0.3,0.6) or (50,150) Myr (early, late) 



for c = (10c, c/10). Evidently, the conse quences of decreas- 
ing t he speed of light (as done by e.g. lAubert fc Tevssierl 
2008) are more dramatic than the consequences of boosting 
it, although neither approximation changes the qualitative 
topology of reionization. 

iMcQuinn et al.l i|2007ft investigated a variety of emissiv- 
ity fields, including cases that were much more biased than 
ours (their models S3 and S4) . Although they did not specif- 
ically address the onset of OI reionization, we may speculate 
on whether their results would have agreed with ours. Their 
fiducial model SI corresponds to a less biased field than 
ours partly because it includes the contribution of all haloes 
down to 10 s Mq, and partly because they assume a constant 
mass-to-light ratio. Their model S3 includes a stronger scal- 
ing between halo mass and luminosity than we find, but it 
still includes the contribution of the low-mass (< 10 8 5 M Q ) 
haloes. Their model S4 assumes that haloes with masses be- 
low 4 x 1O 1O M0 do not form stars, hence it corresponds to 
an extremely biased emissivity field. However, because they 
tune the emissivities to match the volume- averaged emis- 
sivity in their SI model, reionization is already well under- 
way before the density contrasts are large enough for the 
filaments to self-shield efficiently. Hence, of their models, 
the S3 would have been the most similar to ours. Unfor- 
tunately, their spatial resolution is slightly lower than ours 
(367/i _1 kpc versus 250/i 1 kpc), which artificially suppresses 
the density contrast and hence the onset of the OI phase 
(Figure [6} . 

It is possible that increasing our mass resolution would 
boost the star formation rates in low mass haloes at early 
times, leading to a less biased emissivity and more ex- 
tended reionization history. However, even in this case 
our results would remain sensitive to the unknown ioniz- 
ing es cape fraction, which may decrease w ith decreasing 
mass dGnedin et al.l 120081; IWise fc Cenl 120091 ) or increasing 
gas fraction (|Oev et al.ll2007l ). Our finding that the topology 
of reionization depends sensitively on assumptions regarding 
baryonic physics in low-mass haloes emphasizes the need for 
improved observational constraints on the faint end of the lu- 
minosity function as well as on the escape fraction of ionizing 
photons. It also emphasizes the need for high-resolution ra- 
diative hydrodynamic simulatio ns of galaxy forma tion dur- 
ing the reionization epoch (e.g. IWise fc Cenll200"9T ) in order 
to tune the assumptions that go into larger-scale simula- 
tions. 

Our limited cosmological volume may introduce some 
uncertainty into our results. It has been shown that the 
lack of long- wavelength density fluctuations sampled by our 
16/i _1 Mpc volume delays overlap by Az = 0.1 and in- 
troduces a n uncertainty of Az = 1 owing to cosm ic vari- 
ance (e.g.. iBarkana fc Loeblliool llliev et alj|2006al ). While 
it has not been shown that these effects qualitatively change 
the topology of reionization, we have explored this possi- 
bility by computing reionization using snapshots extracted 
from a 32/i _1 Mpc cosmological hydrodynamic simulation. 
The baryonic physics treatments in this simulation are the 
same as in our fiducial volume, as is the assumed ioniz- 
ing escape fraction. Its spatial resolution is 8x lower, with 
the result that star formation is artificially suppressed in 
haloes below the 64-particle limit of 1.5 x 1O 9 M0. We use a 
grid of 64 3 cells, making Axn identical to the tests in Fig- 
ure [9] The dot-dashed magenta curve in Figure [9] shows the 
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Figure 9. Systematic effects on the median reionization redshift as a function of overdensity. The horizontal segments on the right 
indicate the corresponding overlap redshifts. Short-dashed blue, solid black, and long-dashed red curves show the results assuming three 
different values c' for the speed of light. The dotted green curve shows the results when the emissivity is derived from the star formation 
rates. The dot-dashed magenta curve uses a cosmological volume eight times as large as our fiducial volume. Because these computations 
used somewhat less accurate Eddington tensors in order to save time, we also include the result at our fiducial accuracy (light grey 
curves). Changing the speed of light or the volume does not qualitatively change the topology of reionization, but it does impact the 
overlap redshift. Changing how we derive the emissivity field does not affect our results either. 



resulting reionization topology. Comparing it to the solid 
black curve, we find that the topology is qualitatively un- 
changed. In detail, overdense regions reionize later owing 
to the delay of star formatio n at lower mass resolution (e.g., 
ISpringel fc Hernquistl2003bl ) . This in turn delays the overlap 
epoch by roughly Az = 0.2. Ideally, we would prefer run this 
test on outputs from an even larger hydrodynamic simula- 
tion that used the same mass resolution. Unfortunately, this 
is not possible with our current computing resources. Nev- 
ertheless, the qualitative agreement between our 16ft _1 Mpc 
and 32/i _1 Mpc volumes suggests both that our topology 
does not owe to cosmic variance, and that using an even 
larger volume would not change our results. It would be 
interesting to compute the reionization of a 16/i~ 1 Mpc or 
32/i Mpc volume using a testbed that is known to yield a 
purely inside-out reionization topology for larger volumes; 
however, this is beyond the scope of our present work. 

An additional source of uncertainty in our results is the 
unknown impact of structures on size scales below our spa- 
tial resolution such as virialized haloes with masses in the 
range 10 5 7 Mq (i.e., minihaloes). Given that the timescale 
for ionizing photons to bore holes through the IGM into 
the voids is comparable to the timescale over which the 
Universe emits enough photons to ionize each baryon once, 
minihaloes — which are concentrated in filaments — could tilt 
the competition further in favor of an IOM topology through 
their ability to d elay the completion of reionization by up 



to Az = 2 (e.g ., Barkana fc Loebll2002l ; ICiardi et al.ll2006l ; 



McQ uinn et all 120071 ). Moreover, analytic and numerical 
studies now agree that minihaloes change the topology of 



reionization, decreasing the number of large H II regions and 
increasing the number of small ones (|McQuinn et al.l 120071 : 
iFurlanetto fc Ohll2005l ) at late times. Although we do not 
actually resolve minihaloes, this finding is completely con- 
sistent with the trend that we found in Figure [7] whereby 
the mean free path to ionizing photons in the latter stages 
of reionization is shorter at higher spatial resolution owing 
to the impact of overdense "holdouts". This, however, is in 
the late stages of reionization (xhii.v > 0.5); what about 
during the early stages? Minihaloes cannot substantially af- 
fect the topology of reionization until the mean free path 
to absorption by the IGM is comparable to mean free path 
to pass within a virial radius of a minihalo, which, for a 
Press-Schechter mass function and our reionization history, 
happens after z — 8. By contrast, Figure [3] indicates that 
the strength of the ionizing background in voids is already 
comparable to the volume average by z — 9. We conclude 
that minihaloes are not likely to prevent the IOM topol- 
ogy because they are not numerous enough to obstruct the 
photons' paths into the voids. 



6 SUMMARY 

We have used our new moment-based radiative transfer tech- 
nique to study how cosmological reionization proceeds on 
precomputed grids of density and ionizing emissivity span- 
ning a 16 3 h~ 3 Mpc 3 volume. We extract these grids from 
a cosmological hydrodynamic simulation that resolves all 
star formation in dark matter haloes more massive than 
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2 x 10 & M Q . This mass threshold accounts for most of the 
relevant star formation prior to the onset of reionization and 
all of it once reionization is well underway since the nascent 
ionizing background suppresses star formation in haloes less 
massive than ~ 1O 9 M0. 

We find that reionization proceeds rapidly from over- 
dense regions into voids owing to the porous nature of the 
IGM, with filamentary structures ionizing last, which we 
call an inside-outside-middle (IOM) topology. This is con- 
sistent with what we obtained previou sly using less sophis- 
ticated methods l|Finlator et al.l 12009), In our models this 
IOM topology arises because filaments cannot self-ionize, 
since haloes less massive than ~a fewxlO 8 M0 do not form 
stars efficiently. While this cutoff may owe to resolution lim- 
itations in our simulation prior to the onset of reionization, 
it is an expected consequence of Jeans mass suppression in 
reionized regions and the primordial cooling floor at 10 4 K. 

The mean density of ionized regions xm/xv drops from 
high values at early times to the cosmological mean den- 
sity following the epoch of overlap, as expected given that 
sources lie predominantly in overdense regions. At higher 
spatial resolution, the tendency for photons to "leak" di- 
rectly into voids suppresses xm/xv more rapidly so that, 
at our highest resolution, the mean density of ionized re- 
gions drops below the cosmological mean density once the 
volume-averaged ionized hydrogen fraction xh ii,v surpasses 
50%. This leaking effect also manifests as a tendency for 
the mean free path for ionizing photons Amfp to increase at 
higher spatial resolution for xnu,v < 10%. At later stages, 
Amfp decreases at higher resolution because the dense con- 
densations that dominate the IGM opacity at late times are 
better resolved. 

Our reionization history underpredicts the integrated 
optical depth to Thomson scattering and overpredicts the 
strength of the ionizing background at z — 6. The low op- 
tical depth could indicate that reionization begins too late 
and that / osc should be increased, while the high ionizing 
background suggests that the IGM opacity at z = 6 is too 
low. Increasing / csc to unity leads to better agreement with 
the observed optical depth while exacerbating the discrep- 
ancy with the observed ionizing background at z = 6. An 
interesting alternative would be to have / osc = 1 for small 
haloes while leaving it at ^ 10% for larger haloes. This would 
increase the electron fraction at early times, when the ion- 
izing background is dominated by low-mass haloes, while 
leaving it unchanged at late times, when the haloes below 
5 x 10 s Mq no longer form stars. This would be analogous to 
models considering Population III star formation and self- 
regulation of star formation in low-mass haloes, which have 
been shown to yield go od agreement with WMAP observa- 
tions (jlliev et alj|2007h . 

The qualitative topology of reionization is robust to 
our choice of spatial resolution, cosmological volume, and 
ionizing escape fraction even though each of these affects 
the overlap redshift. Increasing the spatial resolution delays 
reionization in overdense regions owing to more effective I- 
front trapping while advancing the reionization of voids ow- 
ing to more efficient "tunneling" of photons through small- 
scale soft spots in the IGM. The latter effect also causes over- 
lap to occur sooner. Doubling the length of our simulation 
volume from 16 to 32/i _1 Mpc delays overlap by Az ~ 0.1 
because the delayed onset of star formation at lower mass 



resolution wins over the tendency for overlap to occur ear- 
lier in larger periodic volumes. The topology, however, re- 
mains essentially unchanged. Increasing the ionizing escape 
fraction to unity causes overlap to occur sooner by Az « 2. 
Additionally, although the filaments are still the last regions 
to ionize in this scenario, they do so somewhat sooner with 
respect to the voids owing to the decreased density contrast 
at higher redshift. This suggests that enhancing the ioniz- 
ing emissivity at early times by, for example, increasing the 
ionizing escape fraction from low-mass halos would lead to 
a more IO-like topology. 

We have also explored how approximations regarding 
the speed of light and the accuracy of the Eddington ten- 
sors impact reionization. Broadly, we find that, while these 
approximations do not qualitatively change the topology of 
reionization, they do affect the overlap redshift. We find 
that (increasing/decreasing) the speed of light by a factor 
of 10 (advances/delays) overlap by (50/150) Myr or Az = 
(0.3/0.6). At a resolution of Axr = 250/i _1 kpc, using Ed- 
dington tensors derived in the optically thin limit delays 
reionization by 40 Myr or Az — 0.2, primarily by delaying 
the reionization of filaments. These results emphasize the 
need to avoid such physical approximations when an accu- 
rate calculation of reionization is desired, and also highlight 
an important feature of our accurate RT code that enables 
us to directly check the impact of such approximations. 

Simulations of reionization are in their infancy. Ideally, 
one would like a single model to span first stars calculations 
on sub-kpc scales within cosmological volumes that encom- 
pass the largest ionization bubbles at overlap. Our current 
codes cannot do so, but future improvements such as an 
adaptive RT mesh and its incorporation into the hydrody- 
namical evolution will bring us closer to this goal. Mean- 
while, we can still gain key insights into the topology of 
reionization, and better understand how one must compute 
reionization accurately. This will set the framework for un- 
derstanding the rapidly advancing observations of this final 
cosmic frontier. 
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